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Assessment of an Explicit Algebraic Reynolds 
Stress Model 


This study assesses an explicit algebraic Reynolds stress turbulence 
model in the in the three-dimensional Reynolds averaged Navier- 
Stokes (RANS) solver, ISAAC (Integrated Solution Algorithm for 
Arbitrary Configurations). Additionally, It compares solutions for two 
select configurations between ISAAC and the RANS solver PAB3D. 
This study compares with either direct numerical simulation data, 
experimental data, or empirical models for several different geometries 
with compressible, separated, and high Reynolds number flows. In 
general, the turbulence model matched data or followed experimental 
trends well, and for the selected configurations, the computational 
results of ISAAC closely matched those of PAB3D using the same 
turbulence model. 


1 Introduction 

The capabilities of application-level three-dimensional Reynolds 
averaged Navier-Stokes (RANS) solvers have increased due to the 
rapid growth of the speed and the size of computational resources 
concurrent with advances in turbulence modeling. Requirements of 
more completely resolving wake-boundary layer flows, highly curved 
flows and jet shear flow physics, for example, have provided impetus 
for the implementation of more, higher order, turbulence models 
into RANS solvers. Two equation turbulence transport models, 
such as k-e or k-w, though a step beyond the previous 0-, 1/2- , 
and 1-equation models, still make significant compromises in the 
physical modeling in the flow calculations. 

Many current three-dimensional Navier-Stokes methods utilize the 
higher order turbulence models such as full Reynolds stress model of 
Launder et a/; 1 the nonlinear model by Shih, Zhu, and Lurnley; 2 and 
the explicit algebraic Reynolds stress model (EASM) of Gatski and 
Speziale. 3 These models have had only a minor additional impact to 
the overall computational effort and, for the most part, have been 
an improvement over linear eddy viscosity turbulence models. 

It is important to be able to evaluate the turbulence models within 
a numerical code whose numerical characteristics are well-known. 
In an attempt to quantify potential numerical biases present in the 
code, a code-to-code comparisons is performed. The challenge is to 
be sure that the turbulence model is similarly implemented in each 
code. This level of effort is time consuming but is becoming more 
essential in the current technical environment where an increasing 
number of validation studies are undertaken. 


1 B. Launder, G. Reece, and W. Rodi. 
Progress in the development of a 
Reynolds stress turbulence closure. 

J. Fluid Mech., 68:537-566, 1975. 

2 T.-H. Shih, J. Zhu, and J. Lumley. A 
New Reynolds Stress Algebraic Model. 
NASA TM-166644, 1994. 

3 T. Gatski and C. Speziale. On Ex- 
plicit Algebraic Stress Models for Com- 
plex Turbulent Flows. J. Fluid Mech., 
254:59-78, 1993 
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The purpose of this paper is to describe and assess the imple- 
mentation of a particular explicit algebraic Reynolds stress model 
(EASM) in the RANS solver ISAAC (Integrated Solution Algorithm 
for Arbitrary Configurations). Earlier results using this model in 
the solver PAB3D are published in papers by Hamid, 4 and Carlson. 5 

The first two sections briefly describe the ISAAC code and 
details of the turbulence model equations used herein. The next 
section contains the validation cases of the EASM simulating four 
different flows and comparing them with the DNS (Direct Numerical 
Simulation), empirical or experimental data where appropriate. 
The last section are N-version tests (code-to-code comparisons) of 
numerical solutions of ISAAC compared with PAB3D, using two 
of the four test cases. The nomenclature is found in Appendix A 
followed by a discussion on turbulence transition indicators used for 
this study. 

2 Numerical approach 

Numerical algorithm 

ISAAC 6 is a multiblock, upwind, finite volume code that solves 
the three-dimensional Favre averaged compressible Navier-Stokes 
equations with an upwind finite volume formulation. The inviscid 
terms are upwind-biased spatial differenced and the viscous terms 
are centrally differenced. This study also uses the flux-difference 
splitting procedure of Roe 7 with MUSCL 8 (Monotone Upstream- 
center Scheme for Conservative Laws) differencing and it advances 
solutions in time with an approximately factored three-factor scheme. 
The limiter used in the MUSCL scheme was Venkatakrishnan’s 
limiter, 9 unless otherwise noted. 

k-e equations transport equations 

The transport equations for the turbulent kinetic energy, k = 
tju'-u 7 , where u 7 is the fluctuating velocity vector, and the turbulent 
dissipation rate, e, are 10 

(pk) =-pP k -pe + V k (1) 

§- t (pe) = C - C £2 f 2 ^ + V £ (2) 

with 

C ei = 1.44, C £2 = 1.92, a k = 1.0, u e = 1.3 (3) 

and where D/Dt = d / dt + Uid / dxi. The equation for the production 
of turbulent kinetic energy, V k , uses the Reynolds stresses and the 


4 K. Abdol-Hamid. Implementation 
of Algebraic Stress Models in General 
3D Navier Stokes Method (PAB3D). 
NASA CR-4702, December 1995. 

5 J.-R. Carlson. Prediction of High 
Reynolds Number Flow Using Alge- 
braic Reynolds Stress Turbulence Mod- 
els, Part 1: Incompressible Flat Plate. 
J. Propulsion and Power , 13:610-619, 
1997; and J.-R. Carlson. Prediction 
of High Reynolds Number Flow Us- 
ing Algebraic Reynolds Stress Turbu- 
lence Models, Part 2: Transonic Shock- 
Separated Afterbody. J. Propulsion 
and Power , 13:620-628, 1997. 


6 J.H. Morrison. A Compressible 

Navier-Stokes Solver with Two- 
Equation and Reynolds Stress Turbu- 
lence Closure Models. NASA CR-4440, 
May 1992. 


7 P. L. Roe. Approximate Riemann 
Solvers, Parameter Vectors, and Dif- 
ference Schemes. J. Comp. Phys., 
43:357-372, 1981 

8 B. van Leer. Towards the Ultimate 
Conservative Difference Schemes V. 

A second Order Sequel to Godunov’s 
Method. J. Comp. Phys., 32:101-136, 
1979 

9 V. Venkatakrishnan. Convergence to 
Steady State Solutions of the Euler 
Equations on Unstructured Grids with 
Limiters. J. Comp. Phys., 118:120-130, 
1995 

10 Several variations of the dissipation 
rate transport equation exist in the lit- 
erature. Many of the variations involve 
the destruction of dissipation term and 
with the singularity of the term as k 
approaches zero at the wall. The cur- 
rent model uses the model proposed 
by Hanjalic and Launder. (K. Han- 
jalic and B. Launder. Contribution 
towards a Reynolds-stress closure for 
low-Reynolds-number turbulence. J. 
Fluid Mech., 74:593-610, 1976.) 
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mean flow gradients, 


V k = -Ta 


dui 


”d Xj 


( 4 ) 


The term T> k derives its form from the turbulent convection and 
pressure transport terms of the Reynolds stress equations. Often, as 
is the case presently, the gradient-diffusion model of Daly -Harlow 11 
is used resulting in 


n = v- 



( 5 ) 


The dissipation rate transport equation uses a fully modeled term, 
T> e , for the viscous diffusion and turbulent transport of e. 


V £ = V- 



( 6 ) 


The eddy viscosity, jit, in these two transport equations is the high 
Reynolds number form with the equilibrium model coefficient , 

l-k = C ,^; Cfio = 0.09 (7) 

The functions in the dissipation rate equation are 

f, = 1, f 2 = 1 - 0.3e- Rt2 , R t = — (8) 

fie 

and the singularity in the destruction of dissipation term is removed 
through the term e 


£ 


£ — 2v 


VVk 


( 9 ) 


The boundary conditions enforced on the wall for the turbulence 
equations are 


k\v — 0 , 


= 2u 




( 10 ) 


Freestream levels for k and e are derived from the ISAAC user 
input variables TKEINF 12 ’ 13 and RMUTNF = /x t //u. TKEINF is then 
nondimensionalized by the square of the speed of sound and used 
as the freestream value for k. RMUTNF determines the freestream 
value for e via the expression for the turbulent eddy viscosity. 14 The 
nondimensional forms designated by the ' are then, 


*4 = TKEINF x Ml = (11) 

£./ _ ^'/.iPoo^'oo f ^ \ Poo Kc f ^ A (V? ) 

00 “ RMUTNF \MocJ ~ (im/v^oo V^oo/ 


11 B. Daly and F. Harlow. Transport 
Equations of Turbulence. Phys. Fluids 
B, 11:2634-2649, 1970. 


12 It is interesting to note that if the 
turbulence intensity is related to the 
ratio of the turbulent kinetic energ y to 
the mean flow kinetic energy, I ~ y/k/F 
(where k = ^u' ■ u', E = -lU-U), then, 
assuming in the freestream that U^q = 
U • U, an expression for the freestream 
turbulent kinetic energy can be written 
as fcoo ~ 


13 More often, and as is the case for 
this study, the streamwise turbulence 
intensity is defined using the stream- 
wise mean and fluctuating velocities 
as I ~ yj u'u' /U. For isotropic flow 
k reduces to ^u'u' , and again solving 
for k and expressing the velocity in 
terms of the Mach number and the 
freestream speed of sound, results in 
koo = 1 The difference 
comes from the definition of I and the 
method of introducing k back into the 
equations. Typical levels for loo are ap- 
proximately 10 — 2 for grid turbulence 
and 10“ 1 for shear flows. The input 
variable TKEINF is equivalent to 


14 fit — R = Reynolds number 

per unit length, and Moo = freestream 
Mach number. 
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Algebraic Reynolds stress model 


One can express the Reynolds stress transport equations in terms of 
the stress anisotropy tensor b (= r /2k — 1/3), 

, Db Dt t Dk 

2k = 

Dt Dt k Dt 

= (v~ T -v^ + n + {v- T -v k ) (is) 

where T> represents the combined terms of viscous diffusion and 
turbulent transport of the Reynolds stresses; T> k = {X>}/2 similarly 
denotes the viscous diffusion and turbulent transport of the kinetic 
energy equations. The braces, { }, denote a tensor inner product. 
II represents the stress redistribution tensor, here modeled by the 
pressure-strain relation developed by Speziale, Sarkar and Gatski 
(SSG). 15 

n = -C x e b + C 2 kS + C 3 k(hS + Sb - ^{bS}l) 

-G^bW-Wb) (14) 

with the symmetric and antisymmetric strain rate tensors, respec- 
tively, 


g _ i /dm 

2 \dxj 


d u. 


W = 


1 ( dm du n 


dxi J 2 \ dxj dxi 

The coefficients of the SSG model are denoted C x — C 4 : 

C, = C{— + Cl, Cl = 3.4, C{ = 1.8, 

C 2 = 0.36, C 3 = 1.25, C A = 0.4 


(15) 


(16) 


Similarly, an expression for the production, "P, in terms of the 
tensors, b, S and W is, 

V = --kS - 2k (bS + Sb) + 2k (bW - Wb) (17) 
3 

Substituting eqs. (17) and (14) into eq. (13) the implicit anisotropy 
transport equation is 

2k 17t h = (' ~ 2Vk ~ Cl£ ^ h c ^ kS 

-(2-c 3 )k(bs + sb + ±-n) 

+ (2-C 4 )l(bW-Wb) + (D-^) (18) 

One can form various models for T > , both differential and algebraic, 
( e.g . the gradient-diffusion model of Daly-Harlow 16 or algebraic 


15 C. Speziale, S. Sarkar, and T. Gatski. 
Modeling the pressure-strain correla- 
tion of turbulence: an invariant dynam- 
ical systems approach. J. Fluid Mech., 
227:245-272, 1991. 


16 B. Daly and F. Harlow. Transport 
Equations of Turbulence. Phys. Fluids 
B, 11:2634-2649, 1970. 
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forms that model T> = /(t, b, T> k , k) 17 ). The present work uses the 
algebraic model for the turbulent transport and viscous diffusion 
by Gatski & Speziale and the solution of the explicit equation 
by Girimaji. 18 The standard algebraic model for 27 , 19 eq.(19), is 
substituted into eq.(18) 

V = jU, (19) 

so that the trailing term is now identically zero. Next, one must 
derive the expression for the production of turbulent kinetic energy 
in terms of the anisotropy and symmetric strain rate tensor from 
one-half the trace of eq. (17), 

n = \{v} = — 2&{bS} (20) 


Substituting in the expression for V k . eq. (20), and assuming 
an equilibrium condition of Dh/Dt = 0, one writes the implicit 
equation for the anisotropy tensor, 


0 = -(c i °-2(2 + C 1 1 ){bS}T)^-a 1 S 

-a 3 (bS + Sb-^{bS}l) |« 2 |bW-Wbj (21) 


a i 
a 4 


Ca), a 2 — -(2 — C 4 ), a 3 — -(2 — C 3 ) 


n , : 

7o— +7i 
£ 


T = 


k 


( 22 ) 

(23) 


One can represent the anisotropy tensor by an expansion of a set of 
basis tensors with scalar coefficients a n as, 


with 


b = «- T(n) 

n = 1 


rp(n) 


S , n = 1 

< s W - WS , n = 2 

^S 2 -i{S 2 }I , n = 3 


(24) 


(25) 


A basis set of three is able to fully describe two-dimensional 
flows. 20 Substitution of eq. (24) into eq. (21) results in a cubic 
equation in scalar coefficient aq, i.e. , 

^) 3 +p(^-) 2 +,(^-)+r = 0 (26) 


17 T. Gatski and C. Speziale. On 
Explicit Algebraic Stress Models for 
Complex Turbulent Flows. J. Fluid 
Mech., 254:59-78, 1993; J.-R. Carl- 
son, N. Duquesne, C. Rumsey, and 
T. Gatski. Computation of Turbulent 
Wake Flows in Variable Pressure Gradi- 
ent. Computers & Fluids , 30:161-187, 
2® 61 -Girimaji. Fully Explicit and Self- 
Consistent Algebraic Reynolds Stress 
Model. Theoret. Comp. Fluid Dynam- 
ics , 8:387-402, 1996 


19 Most explicit algebraic Reynolds 
stress models make this assumption 
for the viscous diffusion and turbulent 
transport terms, T> = , of the 

stress anisotropy transport equations. 
A good discussion on algebraic stress 
models and explicit representations can 
be found in T. Gatski and T. Jongen. 
Nonlinear eddy viscosity and algebraic 
stress models for solving complex tur- 
bulent flows. Progress in Aerospace 
Sciences , 36(8):655-682, 2000. 


20 An evaluation of 3- vs 5-term basis 
sets can be found in T. Jongen and 
T. Gatski. General Explicit Algebraic 
Stress Relations and Best Approxi- 
mation for Three-Dimensional Flows. 
IJES, 36:739-763, 1998. 
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where 21 


V = 


q = 


r = 


7i 

»7i7o 

1 

(2i?i7o) 2 

7i q i 

(27t7o) 2 


7 1 - 277*700! - g^a| - 2?/^ 


70 = ^(Cf + 2), 7i = ^(C' 1 °-2) 

r/t = T 2 {S 2 }, rf 2 = T 2 {W 2 } 

There are three roots to eq. (26) and can be written as, 


P 


77. i = — - + 2 cos - 


P 



„ ^ . e 2 vr\ 

77-2 — — 77 + 2 COS ( 77 H — — J 

9 47t\ 



P 



H 3 — — 77 + 2^/ — T7 C0S I 77 + “g" ) 


with 



b = 2p 3 - 9q 2 + 27r 



6 = arccos(cj)), cf> = 


- 6/2 

V / -« 3 /27 


(27) 

(28) 

(29) 

(30) 

(31) 

(32) 

(33) 

(34) 

(35) 

(36) 

(37) 

(38) 


The only physically realizable solution comes from the smallest 
negative root of the cubic equation. An examination of the 
equations when d < 0 and whether 6 is either positive or negative 
reveal that 77-2 will always be the smallest negative root using the 
following logic: 

If d < 0, 77-2 will always be the smallest negative root, therefore 


aq 


C* = - 77^ = -H 2 


T 


(39) 


If d > 0, there are two roots, 


Hi = ~^P + sign + Vd^j - ^ + Vd 


1/3 


+sign - Vd\ 


4 


1/3 


(40) 


21 The first two invariants of S and W 
are usually defined as 771 = {S 2 } and 
772 = {W 2 }, respectively. Some previ- 
ous articles have used the nomenclature 
of rj 2 (= 771) and U (= -772/771). 
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1/3 



where the model coefficient is 

C* = -(^) = -min(K 1 ,'R 2 ) (42) 

Oil = O 2 O 4 CII (4:3) 

= — 2a3a4«i (44) 


One then writes the Reynolds stress relation as, 


r = -kl + 2/cai T (1 ) + 2 ka 2 T 1 2 3 4 ) + 2 ka 3 T® 
3 

2 2C* AT 2 

= -kl - 2C* fcTS - , ^ / L, 

3 M 4 7o C*T 2 + 7l 


02 


SW-WS) -2a 3 (s 2 -^{S 2 }l)) 


(45) 


(46) 


In order to capture the boundary layer log-layer region, one uses the 
expression C* in place of C* in the linear stress term, making the 
final expression for the Reynolds stresses 


2 2C*fcT 2 

r = -kl- 2C*f„fcTS - , ^ _ 

3 M M 4 7 o C*T 2 + 7l 


H 


S W - ws - 


) -2a 3 (s 2 -i{S 2 }l)) (47) 


3 Results and Discussion 

This section presents four configurations concerning the initial 
validation of this EASM in the CFD (computational fluid dynamics) 
code ISAAC. 


1. Channel flow DNS at R r = 590. 22 

2. Zero pressure gradient flat plate flow DNS at Re = 1410. 23 

3. Backward facing step experiment . 24 

4. Transonic flow over an axisymmetric afterbody experiment . 25 

This study will also compare solutions utilizing the ISAAC flow 
solver of the the zero pressure gradient flat plate and the transonic 
flow over the axisymmetric afterbody with solutions from using the 
solver PAB3D. 


22 R. Moser, J. Kim, and N. Mansour. 
DNS of Turbulent Channel Flow up to 
Re r = 590. Phys. Fluids , 11:943-945, 
1999. 


23 P. Spalart. Direct Simulation of 
Turbulent Boundary Layers up to Reg 
= 1410. J. Fluid Mech ., 187:61-98, 
1988. 


24 D. Driver and H. Seegmiller. Fea- 
tures of a Reattaching Turbulent Shear 
Layer in Divergent Channel Flow. 
AIAA J., 23:163-171, 1985. 


25 D. Reubush. Experimental Investiga- 
tion to Validate Use of Cryogenic Tem- 
peratures to Achieve High Reynolds 
Numbers in Boattail Pressure Testing. 
NASA TM X-3396, August 1976. 
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ISAAC validation test cases 
Turbulent channel flow at R r = 590 


DNS data — Computational results from the RANS code will be 
compared with the numerical experiment (DNS) of Moser et al . 26 
for a fully developed channel flow at a Reynolds number (Re T ) of 590. 
The direct numerical method used a Chebychev-tau formulation 
in the wall-normal direction and a Fourier representation in the 
two other directions. The DNS calculation used a third-order 
Runge-Kutta time discretization for the nonlinear terms, and it 
applied periodic boundary conditions in the streamwise and spanwise 
directions. The computational space was dimensioned 384x257x384 
with a streamwise length of 2n5h, ( 5h = channel half-height ). 

Computational mesh for ISAAC — The grid consisted of a single 
block of an H-type mesh topology with the dimensions of 201x201 
with the streamwise direction in i and the normal (viscous) direction 
in j, respectively. Figure 1 shows the one-quarter density grid. The 


farfield 



26 R. Moser, J. Kim, and N. Mansour. 
DNS of Turbulent Channel Flow up to 
Re T = 590. Phys. Fluids, 11:943-945, 
1999. 


(1/4 grid density) 

Figure 1: Channel flow computational grid with boundary conditions. 

inflow boundary was -0.5 units upstream of the channel entrance. 

The channel length was 5 units, the channel full height (i.e., from 
wall to wall) was 0.01 units, and the first grid point away from 
the wall was at 5.0 xlCP 6 units (the first cell center was placed at 
y + ~ 0.3 ) and geometrically stretched away from the wall at a rate 
of approximately 4 percent. The grid spacing at the channel leading 
edge was 0.025 units and stretched downstream from the leading 
edge at a rate of approximately 0.5 percent. 

This study applied a ‘subsonic inflow’ boundary condition to the 
inflow face i = 1, setting the total temperature and total pressure to 
the input freestream conditions; additionally it applied a ‘subsonic 
outflow’ boundary condition to the outflow face i = irnax, with the 
back pressure to the reference static pressure, p/poo, set to 1. Since 
the subsonic inflow boundary cannot immediately abut a no-slip 


wall, the first 40 cells of the upper and lower boundaries had a 
farfield (1-D Riemann invariants) boundary condition imposed. 

This study also set the initial freestream turbulent eddy viscosity 
ratio (RMUTNF) to 10, and turbulence intensity (TKEINF) to 1.x 10 -4 . 
The freestream Mach number is 0.20, and the Reynolds number per 
unit grid length is 5.8xl0 6 in order to achieve a channel Reynolds 
number, R r , based on the local friction velocity, u r , of 590 in the 
channel. The input unit Reynolds number was manually iterated 
until the correct R r developed in the channel. 

Grid density study — This study calculated solutions at four different 
grid density levels. Each coarser grid was generated by removing 
every other grid point from the previous level grid. 

Figure 2(a) shows the variation with the grid density of the 
channel Reynolds number versus local Reynolds number. 



(a) R r vs. Ri ; ■ 
1/4 grid ; — • • - 


R 

X 

, Baseline grid ; — 
1/8 grid 


& 


800 r 

700 : 
6oo : 
500 : 
400 : 

300 


10 


R = 590 


10 

N 


cell 


10 s 


1/2 grid ; ,(b) R T vs. N ce ;; ; O, R^ = 2.0 x 10 7 ; R x = 5.0 x 10 6 

Figure 2: Effect of grid density. 


With the exception of the region of developing flow (R x < 7 x 10 6 ), 
the 1/4, 1/2 and baseline grids all result in very similar R T 
distributions. The channel Reynolds number at two streamwise 
stations, R x = 5.0 x 10 6 and R x = 2.0 x 10' is plotted against 
the total number of computational cells, fig. 2(b). The channel 
Reynolds number at each streamwise station in the region of 
developing flow, R x = 5.0 x 10 6 , appears to have some local 
grid sensitivity changing roughly 1 percent between the 1/2 and 
baseline grids. The channel Reynolds number at the station further 
downstream at R x = 2.0 x 10' varies less than 0.3 percent between 
the 1/2 and baseline grids. 
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(a) Law-of-the-wall 


Figure 3: Effect of grid density at 
1/4 grid (R r = 586) ; <, 1/8 grid 



(b) Turbulent kinetic energy 

: T = 590: , Baseline grid (R r 

R r = 565). . 



(c) Turbulent shear stress 

590) ; □, 1/2 grid (R T = 588) ; >, 


Mean flow, turbulent kinetic energy, and turbulent shear stress 
are plotted with grid density in fig. 3. All of the profiles show good 
grid convergence. Very little shift is observed between solutions 
generated using the 1/4, 1/2 and baseline grids. 



(a) Mean flow (log-linear scale) 


(b) Turbulent kinetic energy 


(c) Turbulent shear stress 


Figure 4: Comparisons with DNS at R r = 590, linear scale. , EASM; O, DNS [Moser, Kim & Mansour]. 


EASM comparisons with DNS — Figures 4-7 show an evaluation of 
the EASM for predicting mean and turbulent flow characteristics of 
the channel flow. In general, there is good agreement between the 
EASM and the DNS data, with the exception of the anisotropy of 
the normal stresses approaching the wall and the channel centerline. 

The normal Reynolds stresses are equal to 2/3 k with linear 
two-equation eddy viscosity turbulence models. The pressure-strain 
relation in the Reynolds stress transport equation (the basis for 
the EASM) can redistribute turbulent kinetic energy between the 
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(a) Mean flow 

Figure 5: Comparisons with DNS 
Mansour]. 



(b) Turbulent kinetic energy 

R t = 590, log-linear scale. 



(c) Turbulent shear stress 

, EASM; O, DNS [Moser, Kim & 


different stress components resulting in anisotropy in the flow as one 
may observe in figs. 6a and 6b. Though in the region approaching 




(a) Linear scale (b) Log-log scale 

Figure 6: Comparisons with DNS at R r = 590, turbulent normal 
stresses. , EASM; O, DNS [Moser, Kim & Mansour]. 


the wall, the stress anisotropy due to the presence of the solid 
boundary in not correctly modeled. The normal Reynolds stresses 
predicted by the CFD begin to depart from the DNS around y + ~ 
100 and, approaching the wall, become isotropic by y + ~ 0.5 (see 
fig. 6b). 

Figure 7 shows a comparison of the predicted turbulent kinetic 
energy budget with DNS. As a result of the damping of the linear 
terms in the expression for the Reynolds stresses (see eq. (47) and 
the location of the maximum and general profile of the production 
of turbulent kinetic energy, V k , the EASM fairly closely matches the 
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Figure 7: Comparisons with DNS at R r = 590, turbulent kinetic energy 
budget. , EASM ; O, DNS [Moser, Kim & Mansour]. 


DNS calculation. As a consequence the mean velocity profile is also 
well matched as seen in fig. 4a. The dissipation rate, e, balances 
the viscous diffusion of k, T> k , at the wall although the wall values 
of these two quantities are elevated. This is likely due to elevated 
near- wall asymptotic value of k (see fig. 5b) . 

The absence of a damping function in the turbulent eddy viscosity 
in the k and e equations contributes to the excessively high level of 
the turbulent transport of k, T&, compared to the DNS. Additionally 
it contributes to the T^ profile peak occurring too close to the wall 
(see EASM at y + « 5, fig. 7). 

Zero pressure gradient flat plate 

DNS data — Spalart 27 used a fully spectral method in space based 
on Fourier series in the directions parallel to the plate and an 
exponential mapping in the normal direction. The time integration 
uses a Runge-Kutta scheme for the transport terms and a Crank- 
Nicolson scheme for the viscous terms. The following section uses 
mean and turbulent flow data at the Reynolds number based on 
momentum thickness of 1410 from this DNS. 

Computational mesh for ISAAC — The 5-unit long, flat plate, single 
block grid had an H-type mesh topology. The computational domain 
included an inflow block extending 2.5 units upstream from the 
leading edge of the 5 unit flat plate. The initial streamwise grid 
spacing at the leading edge of the plate was 5.0 x 10~ 4 units and was 
exponentially stretched from the leading edge to the trailing edge 
at a rate of 4 percent with a total of 161 grid points. The first cell 


27 P. Spalart. Direct Simulation of 
Turbulent Boundary Layers up to Reg 
= 1410. J. Fluid Mech., 187:61-98, 
1988. 
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height of the baseline grid was 1.0 x 10~ 6 units ( first cell center at 
y + « 0.1 ) fixed at both ends of the plate and exponentially stretched 
from the surface to the outer boundary at a rate of 6 percent with a 
total of 201 grid points. The upper boundary was 5 units away. The 
baseline grid had the dimensions 201 x 201. Numerical transition 
to turbulent flow occurred around R x ~ 0.5 x 10 6 or R e « 500, 
which corresponds to a physical distance of approximately 0.004 
units downstream of the plate leading edge. Grid cell counts were 
divisible by eight to allow three levels of mesh sequencing and/or 
the use of a multigrid to assist solution convergence. 

Grid density study — Figure 8 shows grid density effects for the 
zero pressure gradient flat plate. The solution functionals, local 
skin friction and boundary layer displacement thickness, are plotted 
against total cell count in fig. 8. These solution functionals are not 
monotonically grid converged at Re = 1410, circle symbols in fig. 
8a, because the skin friction of the 1/2 density solution is slightly 
higher than both the 1/8 density and baseline grid density solutions. 
The variations between the skin friction at Re = 1410 for these 3 
finest grids are within 1 percent. One sees a similar trend in the 
variation of the boundary layer displacement thickness at the same 
two stations, fig. 8b, where the 1/2 density grid solution produces 
results below those of both the 1/8 density and baseline density 
solutions. This is due to the point of transition from laminar to 
turbulent flow occurring further downstream for the 1/2 density 
grid than any other grid. 

Further downstream of the turbulent transition point, the 
variation of local skin friction at Re = 1.0 x 10 4 , square symbols 
in fig. 8a, did not change for the 3 finest grids. Variations in 
the turbulent transition point can be masked because the skin 
friction correlation ( Cf vs. Re ) use local functional values, that 
is, the local skin friction, c f = 2 t w / pu^, and a Reynolds number 
based on the local momentem thickness. Changes in the transition 
location globally affect the displacement thickness, hence the same 
non-monotonic grid density trend occurs at both Re stations, see 
square symbols in fig. 8b. 

Figures 9 and 10 shows the variation in the solution with the grid 
density for the mean flow, turbulent kinetic energy, Reynolds shear 
stress at (Re = 1410), the skin friction correlation and boundary 
layer thickness. The mean flow solution changed little with increase 
in grid resolution. Slightly more variance is observed in the turbulent 
kinetic energy and the turbulent shear stress at y + ~ 200 and above. 
This could be in part due to the differing states of transition of the 
boundary for each grid resolution solution. 

The survey station of Re = 1410 was only slightly downstream of 
the start of the numerical laminar-to-turbulent transition phenomena 


0.006 r — j 

0.004 - „ :y 

o“ - n „ „ - 

0.002 - 

o L , 1 , 

10 3 10 4 10 5 

N ce „ 

(a) Cf vs. total number of compu- 
tational cells 

0.006 r j 

0.004 - 

to - ° J * * — - 

0.002 - 


O h , , 

10 3 10 4 10 5 

N cell 

(b) 5* vs. total number of compu- 
tational cells 

Figure 8: Effect of grid den- 
sity on local skin friction and 
displacement thickness at fixed 
stations. O, Re = 1410. ; □, 

Re = 1.0 x 10 4 ; , Baseline 

grid reference line. 
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(a) Law-of-the-wall profiles 


(b) Turbulent kinetic energy 


(c) Turbulent shear stress 


Figure 9: Effect of grid density at Rg = 1410. O, Baseline grid ; □, 1/2 grid ; >, 1/4 grid ; <1, 1/8 grid . 



(a) Skin friction correlation 


(b) S* vs R x 


Figure 10: Effect of grid density at Rg = 1410. O, Baseline grid ; □, 1/2 grid 


D>, 1/4 grid ; <1, 1/8 grid . 
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(constrast, the location of Re = 590 in figures 10a and 10b). Some 
distance appears to be required before ‘fully’ turbulent flow is 
achieved; therefore, in this transitional region the profiles of various 
quantities are subject to some distortion. The difference in the state 
of the boundary layer is not as apparent in fig. 10a as it is in fig. 
10b, noting where Rg = 1410 falls in each plot, and that the variance 
in (5* is observed. Appendix B briefly discusses various indicators of 
laminar or turbulent flow. 


EASM comparisons with DNS — Figures 11 through 13 shows 
comparisons between the EASM and DNS for various quantities. 
Here as well, the EASM closely matched the DNS for the mean 
flow, turbulent kinetic energy, and the Reynolds shear stress. The 
general characteristics of the normal Reynolds stresses are captured 
similar to the channel flow, again with the exception of missing the 
nearwall anisotropy. 




+ 

y 



+ 

y 

(c) Turbulent shear stress (linear scale) 

= 1410. , EASM ; O, DNS 


(a) Mean flow (log-linear scale) (b) Turbulent kinetic energy (linear scale) 

Figure 11: A comparison of the mean velocity profiles with DNS at Rg 
[Spalart], . 


The skin friction correlation, fig. 14a, has a distinct laminar 
characteristic followed by an apparently abrupt transition to 
turbulent like flow starting around Rg = 400. The growth of the 
displacement thickness, <5*, displays a more gradual transition from 
the laminar to turbulent like flow spanning the Reynolds number 
range, R x ss 4xl0 5 to roughly 2.5 xlO 6 , reference fig. 14b. 28 The 
trend and level of both the skin friction correlation and the growth 
of displacement thickness duplicate the classical theory well. 


28 &* = l.72$x/y/R x (Blasius) ; 

5* = 0.046xR^ 1 ^ 5 (turbulent); (F. M. 
White. Viscous Fluid Flow. McGraw- 
Hill, Inc., 1974) 


Backward facing step 

Experimental data — Driver et al . 29 obtained the data for this section 
from an experimental investigation that collected flow survey data 
for a rectangular channel with a backward facing step. The channel 


29 D. Driver and H. Seegmiller. Fea- 
tures of a Reattaching Turbulent Shear 
Layer in Divergent Channel Flow. 
AIAA J., 23:163-171, 1985 
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(a) Mean flow (log-linear scale) 

Figure 12: A comparison of the 
[Spalart]. . 



(b) Turbulent kinetic energy (log-log scale) 

n velocity profiles with DNS at Re 



(c) Turbulent shear stress (log-log scale) 

1410. , EASM ; O, DNS 




(a) Turbulent normal stresses (linear scale) (b) Turbulent normal stresses (log-log scale) 


Figure 13: A comparison of the mean velocity profiles with DNS at 
Re = 1410. , EASM ; O, DNS [Spalart]. . 
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(a) Skin friction correlation (b) vs R x 


Figure 14: Comparisons with classical theory. O, EASM ; , Blasius 

; , turbulent. 


consisted of a 1-m long, 10.1-cm high by 15.1-cm wide channel, 
followed by a 1.27-cm high step down. The top wall of the channel 
varied in angle from -2 degrees (inward) to 10 degrees (outward), 
creating either an adverse or positive pressure gradient downstream 
of the step. This study examines top wall angles of 0 and 6 degrees. 

In units of the step height H, the upstream channel was 8 units 
high, followed by a constant channel height of 9 units for the 0 
degree wall configuration. The channel height was 12.2 units for 
the 6 degree wall configuration at the outflow boundary 40 units 
downstream of the step face. In absolute units, the reference 
conditions were U re f = 44.2 m/s, M re f = 0.128, 5gi = 1.9 cm, and 
Rg = 5000 obtained from data at station x/H = -4. 

Computational mesh for ISAAC — The computational domain 
consisted of two blocks; the upstream block modeled the region 
upstream of the backward facing step; and the downstream block 
included the step and region downstream of the step. The upstream 
and downstream blocks had the dimensions 121x113, and 201x225 
respectively, with a 1:1 block interface. The cell height of the first 
cell was 1.0 x 10 -6 m and was stretched away from the wall at the 
rate of approximately 15 percent until reaching the centerline of the 
channel. The streamwise grid spacing at the step wall was 5.0 x 10~ 4 
nr and was stretched away from that point at approximately 5 
percent in the upstream and downstream directions. 

The location of the inflow boundary was placed 1.0 m upstream 
of the step with subsonic conditions applied along that boundary. 
Turbulent flow developed “naturally” along the 1.0 m wall leading 
up to the location of the step. An alternative would be to define 
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inflow boundary height = 8H 
outflow boundary height = 9H 


x=-80 H 


wall x=0 


x=40H 



backstep (height H) 


subsonic 

outflow 


(1/4 grid density) 


Figure 15: Backward facing step grid . 


a set of fixed conditions at the reference station x/H = -4 derived 
from the specified experimental quantities. The backpressure along 

Table 1: Flow Parameters for the Backward Facing Step 


Case 

(p/Poo) outflow 

M re f 

$bl,ref (cm) 

^6, ref 

Reference data 


0.128 

1.9 

5000 

EASM-ct = 0 deg. 

1.00126 

0.131 

1.70 

6160 

EASM-ct = 6 deg. 

1.0065 

0.121 

1.87 

4750 


the downstream outflow boundary was varied to approximate the 
stated conditions at the reference station of x/FI = -4, see table 1. 
Figure 15 shows the entire computational domain and fig. 16 shows 
the detail of the gridding in the vicinity of the step. 

Grid convergence — Figure 17 shows the variation of static pressure 
coefficient and local skin friction with grid density for the a = 0 
configuration. Similar solutions were calculated using the 1/4, 1/2, 
and baseline grids, the 1/8 density grid departing significantly from 
the other three grids. 

Comparisons with experimental data — Figures 18 and 19 show 
comparisons of the experimental data with the ISAAC solutions 
using the EASM for both a = 0 and a = 6 degree configuration. 
The change in the static pressure coefficient with a variation in the 
wall angle is quite well captured though the actual levels are slightly 
elevated. The change in the local skin friction coefficient with the 
wall angle was similarly predicted though not quite matched as well. 
These results are similar to what has been published previously for 
the Rumsey-Gatski EASM for the a = 0 configuration. 30 The earlier 
results predicted the shift in the boundary layer edge velocity with 
the wall angle, fig. 18c, though the absolute levels were low in the 
region downstream of the step probably due to a high back pressure. 
These results also predicted a change in the flow reattachment point 
with the wall angle though the level was shifted upstream of the 



(baseline grid density) 


Figure 16: Detail of step for the 
backward facing step grid. 


30 C. Rumsey and T. Gatski. Summary 
of EASM Turbulence Models in CFL3D 
With Validation Test Cases. NASA 
TM-2003-212431, June 2003. 
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(a) Static pressure coefficient (b) Local skin friction coefficient 

Figure 17: Variation with grid density, backward facing step, a = 0. , Baseline grid ; , 1/2 grid ; 

, 1/4 grid ; , 1/8 grid. 



-5 0 5 10 15 20 25 30 35 
x/H 



-5 0 5 10 15 20 25 30 35 
x/H 



x/H 


(a) Static pressure coefficient (b) Local skin friction coefficient (c) Normalized boundary layer edge velocity 

Figure 18: Variation with wall divergence angle, backward facing step. O, Data, a = 0 ; □, Data, a = 6 
[Driver and Seegmiller] ; , EASM, a = 0 ; , EASM, a = 6. 
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data, fig. 19. The error bar on the experimental data is from the 
published uncertainty in the determination of the point of flow 
reattachnrent. The error bar applied to the CFD data was the local 
grid spacing. The local grid spacing was considered to be the finest 
resolution possible for determining the reattachment point. 

Axisymmetric afterbody 

Test facility — This test case was an axisymmetric geometry that 
was part of a series of models tested in both the Langley 1/3-m 
Pilot Transonic Cryogenic Tunnel, fig. 20, and the Langley 16-Foot 
Transonic Tunnel. The Pilot Tunnel had an octagonal test section 
with slots at the corners of the octagon similar to the Langley 16-Foot 
Transonic Tunnel test section. 31 The test medium for the cryogenic 



Figure 20: Axisymmetric afterbody model in 1/3-meter Pilot Transonic 
Cryogenic Tunnel . 


tunnel was liquid nitrogen cooled air. High Reynolds number 
data were obtained in the Pilot Tunnel through a combination of 
cryogenic freestream temperatures and freestream total pressure 
that were independently controllable. Approximately 5 atm of 
pressure and 100K total temperature produced a unit Reynolds 
number of 260 x 10 6 /m. 

This test case conducted the experiment over a range of tem- 
peratures from approximately 100K to 300K and pressures from 1 
to 5 times the standard atmospheric level. Several variations of 
freestream total temperatures or pressures resulted in the same 
free stream Reynolds number. Surface pressure coefficients and 
nozzle-boattail drag were shown to be similar, regardless of the 
temperature and pressure combinations that created the equivalent 
Reynolds numbers. 

Previously 32 high Reynolds number simulations with a CFD 
solver were obtained through increased total pressure alone rather 
than through a combination of freestream total pressure and using 


12 ^^ 

10 : 

, 8 - 

6 ; i 

4 t 


-4 -2 0 2 4 6 8 10 12 
a 


Figure 19: Variation of reat- 
tachment point with wall angle, 
backward facing step. •, Data 
[Driver and Seegmiller] (Error 
bar indicates stated experimen- 
tal uncertainty) ; O, EASM 
(Error bar indicates local grid 
spacing) . 


31 R. Kilgore, J. Adcock, and J. Ed- 
ward. Flight Simulation Character- 
istics of the Langley High Reynolds 
Number Cryogenic Transonic Tunnel. 
AIAA Paper 74-80, January 1974. 


32 J.-R. Carlson. Prediction of High 
Reynolds Number Flow Using Alge- 
braic Reynolds Stress Turbulence Mod- 
els, Part 2: Transonic Shock-Separated 
Afterbody. J. Propulsion and Power , 
13:620-628, 1997. 
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cryogenic temperatures. The temperature range of Sutherland’s law 
used for calculating bulk viscosity imposes some limits to the input 
flow conditions. ISAAC input parameters are Reynolds number per 
unit grid length, static temperature and Mach number, so similarly, 
the high Reynolds number conditions were obtained through a 
change in density rather than cyrogenic temperatures. 




(a) 1/4 grid. (b) Detail of afterbody region-baseline grid 

Figure 21: Axisymmetric afterbody computational grid 


Axisymmetric Afterbody Model and Conditions — The model vali- 
dation geometry presented herein is one of six models that were 
built for a Reynolds number investigation performed by Reubush. 33 
This model had a characteristic length of 16 in. and the boattail 
geometry was a circular-arc, with a length-to-maximum-diameter 
ratio (fineness ratio) of 0.8 boattail. The nose of the model was 
a 28° cone, 1.7956 in. long, fairing to the cylindrical body via 
a 1.3615-in. -radius circular-arc whose center is 2.125 in. downstream 
of the model nose and 0.8615 in. below the model centerline. 

The circular-arc fairing is tangent at its endpoints to the conical 
nose (1.7956 in. from the nose) and cylindrical body (2.125 in. from 
the nose). The model was sting mounted with the diameter of the 
sting being equal to the model base diameter. The length of the 
constant diameter portion of the sting (6.70 in. measured from the 
nozzle connect station) was such that there should be no effect of 
the sting flare downstream of the nozzle trailing edge on the boattail 
pressure distributions. These numbers will produce an analytically 
consistent surface definition. The model was cast aluminum and 
from conversations with the original test enginer, it was felt that 
the models were possibly not manufactured to tolerances that would 


33 D. Reubush. The Effect of Reynolds 
number on Boattail Drag. AIAA pa- 
per 75-63, January 1975; D. Reubush 
and L. Putnam. An Experimental and 
Analytical Investigation of the Effect 
on Isolated Boattail Drag of Varying 
Reynolds Number up to 130 x 10 6 . 
NASA TND-8210, May 1976; and 
D. Reubush. Experimental Investiga- 
tion to Validate Use of Cryogenic Tem- 
peratures to Achieve High Reynolds 
Numbers in Boattail Pressure Testing. 
NASA TM X-3396, August 1976 . 
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have produced an analytically smooth shape. 

The grid utilizes an H-type mesh topology, see fig. 21. The block 
dimensions are divisible by 4 to facilitate investigating grid density 
effects on the flow solution. The mesh is gridded with the streamwise 
flow direction oriented along the i index and the j index in the 
normal direction. The body is described using 100 cells extending 
from the leading edge of the nose to the nozzle connect station 
and 80 cells from the nozzle connect station to the nozzle-boattail 
trailing edge, as well as, 80 cells extending 80 body radii from the 
body surface to the far field. The boundary layer grid expansion 
rate away from the body is approximately 16 percent. The inflow 
boundary is 40 body radii upstream of the model nose and the 
outflow boundary is 35 radii downstream of the nozzle boattail. 
Solid walls are treated as no-slip adiabatic surfaces. Riemann 
invarinats along characteristics are specified for boundary conditions 
along the inflow and the lateral freestream outer boundaries of the 
flow domain. The extrapolation boundary condition is applied on 
the downstream outflow face. The flux limiter is min-mod and mesh 
sequencing is used to enhance numerical convergence of the solution. 


Table 2: 

Schedule 

of Reynolds 

Number 

and Grid Spacing 

Nor- 

mal to the Wall 



R(x 10 s ) 

R(x 10 6 )/i 

in. hi 

(in.) 

7.0 

0.4375 

6x 

10“ 5 

128.3 

8.0313 

2x 

10" 6 


0.4 
0.2 
0 

pH 

U - 0.2 
- 0.4 
- 0.6 

-10 12 3 -10 12 3 




X/d m X/d m 

(a) R = 7x10® (b) R = 128.5x10 s 

Figure 22: Variation of afterbody pressure coefficient with grid density, M = 0.9, L = 16. O, baseline grid; 
□ , 1/2 grid, >, 1/4 grid . 


The freestream conditions for the axisymmetric CFD cases were 
M = 0.9, Too = 300F, using air at 7 = 1.4. The first cell height 
of the grid of each configuration was different for each freestream 
Reynolds number according to table 2. Input Reynolds numbers of 
4.375 x 10 5 /in. and 8.0313 x 10 6 /in. were used to develop the low 
and high Reynolds number solutions respectively. 
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Figure 22 shows grid sensitivity of the EASM at M = 0.9 at 
the lowest and highest Reynolds number for this test case, (a) R 
= 7 x 10 6 and (b) R = 128 x 10 6 respectively. These sensitivities 
were relatively consistent for the other turbulence models and 
other viscous models investigated. 34 The EASM prediction of 
the afterbody shock strength and pressure recovery was Reynolds 
number dependent, see fig. 23. 


34 J.-R. Carlson. Prediction of High 
Reynolds Number Flow Using Alge- 
braic Reynolds Stress Turbulence Mod- 
els, Part 2: Transonic Shock-Separated 
Afterbody. J. Propulsion and Power , 
13:620-628, 1997. 




x/d 

m 

(b) R = 128.5X10 6 


Figure 23: Comparison of afterbody pressure coefficients with experimental data at M = 0.9, L = 16. 
EASM ; O, 4> = 0; □, <fi = 120; o , <f> = 240 Data [Reubush] . 


N— Version Testing 
Flat plate 

Figures 24 through 30 compare the solutions produced by the 
ISAAC code with those produced by the PAB3D code for the zero 
pressure gradient flat plate. Several flow variables, such as mean 
flow and turbulence quantities, are compared at two stations, Re = 
1410 and Re ~ 10000, followed by comparisons of the development 
of boundary layer displacement thickness, shape factor and the skin 
friction correlation. 

A very good match has been obtained between the two solvers for 
the mean flow at station Re = 1410, see fig. 24, but the turbulent 
kinetic energy and Reynolds shear stress at station Re = 1410 are 
not as close, see figs. 25 , 26 and 27. Further downstream of the 
numerical transients caused by the transition to turbulent flow at 
the station Re ~ 10000, the two codes appear to match more closely, 
see figs. 28 and 29. 

The first boundary layer shape factor, H 12 and the skin friction 
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(a) Mean flow, Log-linear scale 


(b) Mean flow, Detail of log-layer 


Figure 24: Zero pressure gradient flat plate at = 1410. •, ISAAC; O, PAB3D . 




(a) Turbulent kinetic energy. 


(b) Reynolds shear stress. 


Figure 25: Zero pressure gradient flat plate at Rg = 1410, linear scale. •, ISAAC; O, PAB3D . 
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(a) Turbulent kinetic energy. (b) Reynolds shear stress. 

Figure 26: Zero pressure gradient flat plate at Rg = 1410, log-log scale. •, ISAAC; O, PAB3D . 



(a) — ai/T (b) P k /e. 

Figure 27: Zero pressure gradient flat plate at = 1410. •, ISAAC; O, PAB3D . 


25 



(a) Log-linear scale. 


(b) Detail of log-layer. 


Figure 28: Zero pressure gradient flat plate at Rg = 10 4 , mean flow. •, ISAAC; O, PAB3D ; • • • , log-layer. 




(b) Log-log scale 


Figure 29: Zero pressure gradient flat plate at Rg = 10 4 , turbulent kinetic energy. •, ISAAC; O, PAB3D ; 
• • • , log-layer. 
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(a) Momentum thickness with local Reynolds 
number 


Figure 30: Zero pressure gradient 
theory-turbulent. 


3 



io 1 io 2 io 3 io 4 io 5 


Re 

(b) First shape factor 

plate. •, ISAAC; O, PAB3D ; 



(c) Skin friction correlation 


, Linear theory-laminar ; , Linear 


correlation, c f ws R^, are also quite closely matched with only 
minor variances seen in the upstream laminar portion of the flat 
plate, see figs. 30b and 30c. Similarly, the development of the 
the displacement thickness with Reynolds number, fig. 30a, for 
the most part is identical between ISAAC and PAB3D. A slight 
departure occurs in the the laminar-to-tui'bulent transition region, 

3 x 10 5 < Rx ^ 2 x 10 6 , but is relatively local and does not appear 
to affect the downstream match significantly. 

Axisymmetric afterbody 

Using the two codes, fig. 31 compares the static pressure coefficients 
on the afterbody. A slightly higher pressure occurs in the ISAAC 
solution compared with the PAB3D solution for the lower Reynolds 
number case, see fig. 31a just downstream of x/d m = 0.0. The 
high Reynolds number comparison show virtually identical static 
pressure distributions on this portion of the model. 

Initially, the solutions developed using ISAAC had a considerably 
weaker shock because of the type of limiter being applied to 
the inviscid fluxes. Since the PAB3D solutions were developed 
using the min-mod limiter, changing the limiter being used in the 
ISAAC calculations from Venkatakrishnan’s smoothing limiter to 
the min-mod resolved this difference in the solutions between the 
two solvers. 

4 Summary 

The turbulence equations being solved in the three-dimensional 
Reynolds averaged Navier-Stokes (RANS) code PAB3D have 
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Figure 31: Comparison of static pressure coefficients, axisymmetric afterbody, M = 0.9, L = 16. , ISAAC 

; , PAB3D ; O, 4> = 0; □, <p = 120; o, <f> = 240, Data [Reubush]. 


been implemented in the three-dimensional RANS code ISAAC. 
Initial validation of the implementation was accomplished through 
comparisons with data for four test cases: channel flow, zero pressure 
gradient flat plate, backward facing step and an axisymmetric 
transonic afterbody flow. These test cases covered a wide range of 
fluid dynamic situations, including compressible flow, separated flow, 
and high Reynolds number flow. The turbulence equations were 
able to simulate all of these flow situations fairly well. Additionally, 
excellent agreement was achieved in comparisons between solutions 
developed using ISAAC and with solutions developed using the 
solver PAB3D. 
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A Nomenclature 


Symbol 

Definition 

b, b ij 

C* 
C p 

c f 
dm 
H 
H 12 
h, S h 
hi 
L 
M 
N ce Zi 
V 

n 

Reynolds stress anisotropy tensor 
Variable model coefficient, — ai/T 
Pressure coefficient, (p - poo) / qoo 
Local skin friction coefficient, r w / qoo 
Body maximum diameter (1.0 in.) [axisymmetric afterbody] 
Backward facing step height 
First boundary layer shape factor, (8* /0) 

Channel half-height 

Physical height of first computational grid from a wall 
Reference length [axisymmetric afterbody] 

Mach number 

Total number of computational cells 
Production, —Tik(diLj /dx^) — Tjk(dui/dxk) 

Production of turbulent kinetic energy, —\{"P} = Tij ( dui/dxj ) 

p 
q 
R 
Rt 
Rx 
Re 
Rr 
S, S ^ 
T 

Static pressure 
Dynamic pressure 

Reynolds number based on model reference length, (uooL /v) 
Cell turbulent Reynolds number, k 2 /(ue) 

Reynolds number based on distance x, (uooX/V) 

Reynolds number based on momentum thickness, (uooO/is) 
Reynolds number based on wall shear stress, (u T h/z^) 
Symmetric velocity gradient tensor, ( dui/dxj + duj/dxi) /2 
Kinematic time scale, k/e 

U, U i 
U + 

Velocity vector 

Law-of-the-wall coordinate, u/u r 

u T 

w, w i:j 

X 

Local friction velocity, y T w /~p w 

Antisymmetric velocity gradient tensor, ( dui/dxj — duj/dxi) /2 
Streamwise distance 

X, Xi 

Position vector 

Xr 

Point of flow reattachment 

y 

y + 

a 

Vertical distance 

Law-of-the-wall coordinate, y u T / v 
Wall divergence angle (backward facing step) 

a n 

8 

5* 

scalar coefficient 

Boundary-layer thickness, value of y at 0.995 u maa: 
Displacement thickness, J’ 0 5 (l-u/u e )dy 

vl 

V2 

e 

Flow invariant, T 2 {S 2 } 

Flow invariant, T 2 {W 2 } 

Momentum thickness, f Q (u/u e ) (1 - u/u e ) dy 

P 

V 

Laminar viscosity 

Turbulent eddy viscosity, C^ Q ~pk 2 /e 
Kinematic viscosity, fi / p 

p 

T, Tij 
T~w 

<t> 

bl 

ref 

Density 

Reynolds stress tensor 

Shear stress at the wall, ~p w du/dn 

Angular location of pressure orifices, deg. (axisymmetric afterbody) 
Boundary layer 
Reference conditions 

w 

Condition at wall surface 

oo 

Freestream condition 
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B A Comparison of Transition Indicators 

There are several parameters that can indicate the presence of 
turbulent flow, or as the case may be, developing turbulent flow over 
the flat plate. These indicators occur both at a primative level, e.g., 
peak turbulent kinetic energy and dissipation rate level at the wall; 
or at a derived level, e.g., local skin friction, the rate of boundary 
layer growth with local Reynolds number, or boundary layer shape 
factors, to name a few examples. 

Typically, a variation of the local skin friction with the local 
Reynolds number, fig. 32 is plotted, which clearly deliniates the 
change from laminar flow by the rather abrupt discontinuity in the 
curve around R x = 400,000. 

If the Reynolds number based on the momentum thickness of the 
boundary layer is used instead of the local Reynolds number, an 
equally abrupt change in the skin friction correlation from laminar 
flow to turbulent flow levls is observed, fig. 33. At first glance, 
the local skin friction in these two figures appears to attain a fully 
turbulent level almost immediately downstream of the transition 
point. (As noted on the figure, the momentum Reynolds number of 
410 is approximately coincident with the local Reynolds number of 
400,000 and is used as a reference point for the discussion.) 



Figure 34: Displacement thickness with local Reynolds number (detail), 
o, EASM ; , Blasius ; , turbulent (flat plate theory) . 

A slightly different conclusion could be drawn (concerning the 
transition of the flow from laminar to turbulent flow) from a plot of 
the rate of growth in the boundary layer displacement thickness with 
local Reynolds number, fig. 34. The displacement thickness only 
gradually approaches the turbulent level as the Reynolds number 
increases downstream of the transition point marked (R# ~ 410). 35 

Similarly the first shape factor, H 12 = S*/0, shows a clear 
departure from laminar flow, fig. 35, but a rather prolonged region 



R 

X 


Figure 32: Local skin friction 
with local Reynolds number, o, 
EASM. 



Figure 33: Skin friction correla- 
tion. o, EASM ; , Karman- 

Shoenherr empirical curve. 


35 M. Rai and P. Moin. Direct Nu- 
merical Simulation of Transition and 
Turbulence in a Spatially Evolving 
Boundary Layer. Journal of Computa- 
tional Physics , 109:169-192, 1993. 
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Figure 35: First shape factor with momentum thickness Reynolds 
number, o, EASM . 


of developing turbulent flow before the shape factor approaches the 
classically derived level. 

If the “simpler” parameter of peak turbulent kinetic energy is 
plotted against Reynolds number, a less definitive picture appears 
as to the point of transition, fig. 36. For reference, the figures show 
the locations of R x ~ 400000 and R# ~ 410. 

Though still somewhat abrupt at transition, the peak turbulent 
kinetic energy has a region of gradual increase in intensity before 
increasing another two orders of magnitude greater than the level in 
the region of laminar flow, in the current example, starting from a 
value around 10 -6 . 
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Figure 36: Streamwise variation of peak turbulent kinetic energy. 
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